Identifying an Experimental Two-State Hamiltonian to Arbitrary Accuracy 
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Precision control of a quantum system requires accurate determination of the effective system 
Hamiltonian. We develop a method for estimating the Hamiltonian parameters for some unknown 
two-state system and providing uncertainty bounds on these parameters. This method requires only 
one measurement basis and the ability to initialise the system in some arbitrary state which is not 
an eigenstate of the Hamiltonian in question. The scaling of the uncertainty is studied for large 
numbers of measurements and found to be proportional to one on the square-root of the number of 
measurements. 

PACS numbers; 03.65.Wj, 03.67.Lx 



I. INTRODUCTION 

High precision control of quantum systems inevitably 
requires high precision characterisation of the system dy- 
namics. Quantum computers are an example of a device 
requiring especially high precision characterisation, but 
this precision is also required for detailed studies of in- 
teractions in quantum systems. For qubits, this charac- 
terisation is usually performed using state and process 
tomography, where the full density matrix is measured 
for a range of different input states Q, |^ H, Q . An al- 
ternative approach is to directly characterise the Hamil- 
tonian, which then gives the evolution of the system for 
any initial state. This approach is especially useful when 
the system approaches a closed system and therefore its 
dynamics can be treated as purely Hamiltonian. While 
this will not be the case in general, it is an essential 
requirement for constructing a qubit for quantum com- 
puting applications and is approximately true for many 
other systems of interest. 

Tomographic methods typically require measurement 
in several different bases or require the ability to per- 
form rotations around particular axes before the system 
has been completely characterised. In contrast, a general 
procedure developed recently for identifying an arbitrary 
two-state Hamiltonian^ requires measurement in only 
one basis and initialisation in a single known state. The 
requirement for only one measurement basis is especially 
attractive for systems with limited measurement devices, 
for example many solid-state qubits. We build on this 
result by deriving a systematic method to calculate the 
Hamiltonian parameters to any required accuracy from a 
time series of measurement data. 

In previous work the Hamiltonian is assumed to be 
a linear combination of some free evolution Hamiltonian 
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and various control fields, where characterisation requires 
finding both the base Hamiltonian and the dependence 
on the control field. In this work, we take a more prag- 
matic approach to the problem of characterisation of a 
two-state system. We provide a method to answer the 
question 'What measurements must be taken to deter- 
mine the form of a two-state Hamiltonian to a given pre- 
cision?' Assuming that the system evolves under some 
Hamiltonian, which corresponds to a certain 'position on 
a dial' in the laboratory, the parameters for this partic- 
ular Hamiltonian can be determined to some arbitrary 
precision. If the two-state system is to be used as a qubit 
for quantum information processing (QIP) applications, 
the process can then be repeated for some other linearly 
independent Hamiltonian, giving two 'axes' that are suffi- 
cient to construct any arbitrary single-qubit rotation |^ . 
In general, the response of the system to various 'dial 
settings' would be required to construct efficient single- 
qubit gates. To do this, the Hamiltonian parameters and 
their uncertainty would need to be determined for a num- 
ber of points and the response determined. In the case 
of linear response, this becomes completely equivalent to 
the process discussed in reference p but more generally 
will require fitting to an appropriate functional form. 

The basic outline of this method and the relevant equa- 
tions are given in sections imillll and l VII The uncertainty 
in these estimates is then analysed and a series of un- 
certainty relations are given in sections IVIVIII and IVIIII 
which allows the Hamiltonian to be estimated with error 
bounds on all its parameters. Section IIVI covers some 
technical details on the use of the discrete Fourier trans- 
form (DFT) to analyse the time series data and how its 
accuracy can be controlled for this particular application. 
In section IIXI we numerically simulate this method for 
some example Hamiltonians and compare the statistical 
spread of results with the estimated uncertainty, finding 
very good agreement. We also investigate how the accu- 
racy of the Hamiltonian parameters scales as a function 
of the number of measurements. Finally, in section Ixl we 
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FIG. 1: Bloch sphere representation of the state of a qubit 
and its trajectory given an arbitrary Hamiltonian d. If the 
system is not in an eigenstate of the Hamiltonian, the state 
given by the unit vector s precesses around an axis defined by 
d. The components of d are given by the Hamihonian using 
Eq. Q where the |d| gives the angular precession frequency 
around the vector (dx,dy,dz)^ . 



discuss the effect of this scaling on the characterisation 
and operation of single qubit gates for QIP applications. 



II. CHARACTERISING A TWO-STATE 
HAMILTONIAN 

Some insight into the time evolution of an arbitrary su- 
perposition state can be gained by considering the Bloch 
sphere picture for a two-state system. The Hamiltonian 
of an arbitrary two-level system can be written in terms 
of the Pauli matrices, 



H — —— = ^-^{dol + dx<Tx + dy(Ty 



(1) 



where dx ,dy and dz are real constants and do results in an 
unobservable global phase factor which can be ignored. 
If the state of the system is mapped to the Bloch sphere, 
its position in the sphere is the Bloch vector (s) where 
|s| < 1, with a pure state having |s| — 1. The evolution of 
the Bloch vector due to some Hamiltonian {H) will be to 
precess around a unit vector {dx,dy,dz)'^ with angular 
rotation frequency given by |d|. If the system is in an 
eigenstate of the Hamiltonian, the Bloch vector is parallel 
to the axis of rotation and therefore does not precess, as 
expected. This process is illustrated in Fig. ^ 

If the system can be repeatedly initialised in a known 
state and then measured in some basis at progressively 
longer time periods, the trajectory of the Bloch vector 
can be mapped. Assuming this is an idealised projective 
measurement, the sinusoidal variation of the projection 
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FIG. 2: To map z{t) the system must be repeatedly initialised, 
allowed to evolve under the Hamiltonian to be measured H{t) 
and then measured. To map the time evolution of the system, 
the Hamiltonian step is applied for progressively longer time 
intervals {iAt for i = 1,2, ... ,n) where At is the minimum 
controllable time interval and iob ~ nAt is the maximum time 
over which the system is observed. 



onto the measurement axis depends on both the magni- 
tude and direction of the vector d and therefore on the 
parameters in the Hamiltonian. A schematic of this pro- 
cess is shown in Fig. |21 where the minimum controllable 
time interval is given by At and the longest time the sys- 
tem is allowed to evolve is toh giving the total number of 
time points Ng = to\^/ At. This process is then repeated 
Ne times to build up an ensemble average for each time 
point, giving a total of Nt — Netoh/^t — N^Ng mea- 
surements. The true evolution z{t) is then approximated 
by the measured function Zm{t). 

For simplicity, we will use polar coordinates to describe 
both the position of the Bloch vector and the Hamilto- 
nian vector, as illustrated in Fig.^ In these coordinates, 
the Hamiltonian vector is given by d = jdKd^,, dy, dz)^ = 
|d|[sin(6') cos(0), sin(6') sin(0), cos(0)]-^. As the complex 
phase (0) is unobservable in a single two-state system, 
we can set (f> = and therefore align the Hamiltonian 
with the X-axis. If the reference axes are defined based 
on experimental grounds, this can be corrected with a 
trivial rotation. 

If the system is initialised in the state \ip{0)) ~ |0) 
(which corresponds to = </> = or Sq = (0, 0, 1)^), the 
evolution of the z-component of the state vector is 



z{t) = cos(cji) sin2(6') + cos^(6'). 



(2) 



where = |d| (in units such that h — 1), see reference 
or appendix^for an alternative derivation. Determining 
the parameters ui and cos^{9) gives the values of |d|, dx 
and dz. Throughout this discussion, we assume that the 
Hamiltonian is constant in time and that the initial state 
of the system is not an eigenstate of the Hamiltonian, i.e. 
0^0, otherwise the system will not precess. The process 
of characterising the Hamiltonian thus involves measur- 
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ing Zm{t) and analysing it to determine the appropriate 
parameters. 



III. ESTIMATING THE HAMILTONIAN 
PARAMETERS FROM FOURIER COMPONENTS 

Once Zm{t) is determined, the data can be fit- 
ted in the time domain to determine the Hamiltonian 
parameters]^. While this is sufficient for approximate es- 
timates or data containing only a few oscillation periods, 
a more elegant method is to take the Discrete Fourier 
Transform (DFT) of the data z„i{t) and calculate the 
parameters from the Fourier coefficients. This method 
provides both the Hamiltonian parameters and an esti- 
mate of the uncertainty in these values. In order for this 
method to be effective, the Hamiltonian must be con- 
stant in time, or more precisely, the fields controlling the 
Hamiltonian must be stable to higher precision than that 
required for characterisation pll |. 

As z{t) is a pure sinusoid, if Zm{t) consists of an inte- 
ger number of periods of oscillation, its Fourier transform 
~ DFT[zra{t)]) will take on a simple form consist- 
ing of (5-functions at = and v = dzi/p where Vp refers 
to the position of the peaks. Using the definition of the 
inverse discrete Fourier transform DFT^^, z„i{t) can be 
rewritten in terms of the discrete Fourier components for 
the zero [-F'(O)] and peak-frequencies [F{vp)], 

Ns/2 

DFT-i[DFT[z,„(i)]] = J2 i^(j^)e'2'^(''/^=)* 

= i^(0) + F(l/p)e*27r,.,t/Ar, 

= F(0) + 2F(fp) cos(27ri/pV^s) 
z{t). 

In this way, the angle and the angular precession fre- 
quency uj can be determined directly from the Fourier 
spectrum without the need for fitting the data in the 
time domain. 

The effect of a measurement error probability can also 
be included by assuming some probability rj € [0, 1] of 
obtaining the incorrect value from a single measurement. 
This corresponds to a bit-flip error {ax) occurring the in- 
stant before measurement with some probability rj. As- 
suming the Bloch vector always starts at = |0), 
z[t) should reach a maximum of one after each period. 
The measurement error will reduce this maximum, inde- 
pendent of the angle 0,. and can therefore be determined 
directly from the DFT. If we model the effect of this mea- 
surement error as Zm(t) = {l — 2r])z{t), then the following 
equations can be derived, 

V=^—^-Fi^p), (3) 




Time (t) Frequency (v) 



FIG. 3; The left hand plot shows an example of a sampled 
time signal z{t) = [cos(27rt) -I- l]/2 with A^e = 1 (a), 2 (b), 
8 (c) and 500 (d) measurements at each time point, where 
each measurement is a projection onto the (1,-1) axis. The 
corresponding DFT for each signal is shown on the right for 
ly > 0, illustrating the signal-to-noise improvement as more 
measurements are taken at each time point. 



iu = 2nVp/Ns. (5) 

As we can only perform projective measurements onto 
one axis, many measurements are required to accurately 
determine Zm(t), so Nt will typically be quite large. 
Once the time resolution and observation time are cho- 
sen, the measurements for each time point can be re- 
peated until a sufficiently resolved peak is seen in the 
DFT spectrum. An example of this process is shown in 
Fig. 121 for progressively larger numbers of measurements 
at each time point. In this way, the number of measure- 
ments need not be chosen at the start but the experi- 
ment is repeated until a sufficient signal-to-noise ratio is 
obtained. 



IV. DETERMINING THE PRECESSION 
FREQUENCY TO ARBITRARY ACCURACY 

Performing a discrete Fourier transform (DFT) on 
the measurement results immediately places some con- 
straints on the selection of the measurement parameters. 
In order to satisfy the Nyquist sampling criteria, at least 
two sample points for every period of oscillation are re- 
quired to avoid aliasing. This means that some estimate 
for the oscillation period Tprcdict must be known in order 
to guarantee that At < Tpi.cdict/2, though in practice the 
period of oscillation will usually be known approximately 
on theoretical or experimental grounds. 

Conventional DFT theory states that the frequency 
resolution {Aiy) of a DFT signal is the inverse of half 
the total time of the signal, — 2/tob0- This means 
that to resolve the frequency signal we need to observe at 
least two complete oscillation periods, though typically 
many more periods will need to be observed to obtain 
a clearly defined peak in the frequency spectra. For an 
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Time (t) Frequency (v) 



FIG. 4: The left hand plot shows time signal which are trun- 
cated at various time points to produce a net phase difference 
of /S.tf = TT (a), /S.tp = Ti/2 (b) and A</p = (c) between the 
start and end of the signal. The corresponding DFT for each 
signal is shown on the right, where the peak approaches a 
(^-function only for ~ 0. 




t 

p 

FIG. 5: The test function P{tp) used to locate the point at 
which there is zero phase difference between the first and last 
sample point. The amount of the time signal to use in the 
DFT is given by tp and the uncertainty is given by the FWHM 

of P(tp). 



arbitrary signal the frequency resolution of the spectra 
also limits the precision with which one can determine 
the frequency {v ± Av) . The more periods observed the 
more accurately the determined frequency of oscillation. 
Ultimately this will be restricted by the decoherence time 
of the system as decoherence reduces the amplitude of the 
oscillations for long observation times. 

To use Eqs. Q-©, we require that the observation 
time tob is an integer number of periods. To ensure 
this, we need to know the precession frequency to the 
same precision as the time control (Aiy/vp « At/toh)- 
Conversely, if we can guarantee that we have an inte- 
ger number of periods, this will yield the corresponding 
frequency. 

The DFT of a pure sinusoid has some special prop- 
erties in that it only approaches a (5-function when the 
time signal consists of an integer number of periods (there 
is no phase difference between the start and end of the 
signal) 0. If there is some phase difference then the 
DFT has 'leakage' into the other channels, resulting in 
a overall spread of the signal throughout the spectrum. 
This effect is demonstrated in Fig. 0] for example si- 
nusoids having various values for the phase difference 
{Aip = ip(0) — ip{tob)) between the start and end points 
in the time signal. 

Using this information, we can locate the 'minimum- 
phase-point' (MPP) where the difference in phase be- 
tween the start and end of the signal is minimised. This 
amounts to selecting only an integer number of periods 
of the signal. As the period of the signal is not known 
beforehand, the easiest method is to record the data and 
then reprocess it later to ignore some of the data points. 
While this results in throwing away some information, 
the lost data consists of at most one period. 

An effective way of locating the MPP is to compare 
the magnitude of the channel comprising the central fre- 
quency peak F{vp) and its adjoining channels F{vp — 1) 
and F{i'p + 1). When the leakage is minimised, the ra- 
tio of the central channel to its neighbours should be a 
maximum. An example test function which was found to 



perform well with varying levels of noise is 

p., ^ _ m'^p) - Fji^P - 1) - FjUp + 1) 

where once again F{v) is the normalised DFT of the orig- 
inal signal from Zrn{0) to Zm{tp) where toh - T'picdict < 
tp < iob- An example plot of P{tp) is shown in Fig. |S1 
A clear peak is observed at the point where the phase 
of the sinusoid (lys) is an integer multiple of 27r, i.e. 
9?(0) = ^p{tp = 27rm) for some integer m. Once the 
MPP has been determined, the frequency is given by 
LJ — 21:71 1 tp where n is the peak channel number and 
tp is the MPP. The advantage of this method is that the 
MPP can usually be determined to an accuracy of close 
to At/toh and the full- width- half-maximum (FWHM) of 
the function P{tp) gives an estimate for the uncertainty 
of the resulting frequency. 



V. ESTIMATING THE UNCERTAINTY IN THE 
MEASURED QUANTITIES 

For most practical applications, if we wish to estimate 
the parameters of a two-state system, we also need to 
know the uncertainty in those estimates. For the rest of 
the discussion we will use the following notation, x is the 
estimate obtained for some true value x and 5x refers 
to the predicted standard deviation of this estimate. In 
the ideal situation x — 3(5x < x < x + 'i5x, with 99.7% 
confidence. 

As we are determining the parameters of interest from 
the components of the Fourier spectrum, we have a 
straightforward way of calculating the uncertainty from 
the spectral noise. We define the noise spectrum n{v) to 
be the parts of the Fourier spectrum which do not in- 
clude F{±Vp) and F(0). This is a good approximation 
when iob constitutes an integer number of periods and 
therefore F{±Vp) and -F'(O) approach ^-functions. 

In general the noise due to the discrete measurement 
of the system will be a limiting factor in the analysis. 
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though other factors hke noise in the control Hamilto- 
nian will also contribute. The uncertainty in the fre- 
quency will be primarily controlled by the precision in 
the time control of the measurements. Ideally the uncer- 
tainty in the angular frequency measurement should be 
of the same order as the time resolution in the measured 
signal {Slo/lo « Ai/^ob)- In practice a more accurate 
estimate for the uncertainty can be obtained from the 
FWHM of P{tp), as discussed earlier. The uncertainty 
in the angle 69 and the measurement error Srj will be 
primarily limited by the noise level in the Fourier spec- 
trum. Typically, the fractional uncertainty in lo will be 
an order of magnitude smaller than for or 77 as finding uj 
only requires finding the peak location where as the other 
parameters depend on the peak height which is directly 
affected by the spectral noise. 

The uncertainty in the Fourier peaks is given by the 
standard deviation (SD) of the noise spectrum. For 
simplicity we will define 6F = ST)[n(iy)] and 6uj = 
27r/FWHM[P(tp)] so that the resulting uncertainty ap- 
proximates the predicted standard deviation of the pa- 
rameter estimate. Once we have the uncertainty in the 
frequency Suj and the Fourier spectrum SF, using con- 
ventional uncertainty analysis [8| we can derive the ex- 
pressions for the uncertainty in the calculated values. 
Throughout this discussion we will use the standard er- 
ror propagation method [sl where the variance of some 
function w — f{x,y) is given in terms of the variances 
var(a;) and var(7;) and the covariance cov(a;, y) between x 
and y '22i] . In its simplest form, the variance of a function 
can be calculated using 



var(i()) = 



OF 



dX 
+2 



var(a;) 



dF 



dF 



dX 



dF 
W 



dY 
cov{x,y), 



var(?/) 



(7) 



for small variances in the measured parameters. 

Using this approach, the uncertainty in each of the 
calculated quantities in Eq. © and 10} can be estimated 
using the following equations. 



6rj 



rSF, 



(8) 



6A' 



m 



l- 277 -F(0) 



Sr/ 



1-277 



6F' 



and 



(1 - 277)3 

(l-A2)-l/2 



(9) 



(10) 



where A = cos(6'). 

This process results in an estimate and its associated 
uncertainty for the angular frequency uj, rotation axis 9 



and the measurement error 77. A simplistic error analy- 
sis is given here to illustrate the ideas. The use of more 
sophisticated techniques such as maximum likelihood es- 
timation should provide tighter bounds on the estimated 
parameters for a given set of data[3nlQ.illj. 



VI. DETERMINING THE PHASE-ANGLE 
BETWEEN TWO HAMILTONIANS 

The process discussed so far is sufficient to characterise 
a single two-state Hamiltonian, as dy can be arbitrarily 
set to zero. To provide a completely controllable two- 
state system, such as is needed for QIP, a second con- 
trol Hamiltonian is required to implement all possible 
single-qubit rotations. If we consider characterising some 
reference Hamiltonian {Hr), we can use this to define 
the coordinate axes and then consider a second Hamilto- 
nian (-fffe). This provides a second axis to rotate around 
which must also be characterised and the angle (j) be- 
tween these two axes must be determined. To measure 
this azimuthal angle, a different initialisation point must 
be chosen whose Bloch vector is linear independent of 
the original initialisation point. A convenient choice is 
to rotate s around the first axis (d^) until in sits on the 
'equator' defined by 9 = 7r/2. The second Hamiltonian is 
then switched on instead and the qubit precesses around 
dfe. The z-projection of this rotation can then be used 
to determine the angle (j) between the two axes. As d^ 
and dfc have already been completely characterised, the 
entire process can be 'boot-strapped', progressively learn- 
ing more information about the system. Of course, this 
process of measuring different Hamiltonians is equivalent 
to measuring the dependence of a system Hamiltonian 
on the 'settings of a dial' where each Hamiltonian corre- 
sponds to a different value for the input parameters. 

To rotate s onto the equator, starting at Sq we apply 
d,. for a time 



1 

t = — arccos 



cos(2^?j,) 



1 



cos{29r) - 1 



(11) 



which places the system in state Si = [cos(/3), sin(/3), 0]-^ 
where f3 = arctan[— sec{9r)\/ —2 cos{29r)] If we then 
use this as the new initialisation point, the z-component 
of the precession about d^ is given by 



z{t) — C[l — cos{uJkt)] + Dsm{iLjkt), 



(12) 



where C = ^ sm{29k) cos(0 — f3) and D — sm{9k) sm{(j) — 
(3). This procedure can only be apphed if 9r G [f j %-]■ 
If 9r or 9k are not within this range, a more elaborate 
pulsing scheme is required. Once the two axes d^ and d^ 
have been characterised, measuring Eq. (|12|l allows both 
Hamiltonians to be completely reconstructed. 

Using a similar method to the previous section, the 
parameters C and D can be determined from the com- 
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ponents of the Fourier spectrum, 



DFT-i[DFT[z™(i)]] 



E 



F{y)e 



= F{Q) + i^ij(zyp)e''2'^''^*/^» 

~iFi{^Vp)e-"'^'"'-^'/'^' 
= F(0) + 2Fr{vp) cos{2nVpt/N,) 

-2Fi{iyp)sm{2Tnyp/Ns) 
^ zit). 

where Fr and F/ are the real and imaginary parts of the 
Fourier components. As the measurement error of the 
system has already been determined from the measure- 
ments of the other axes, the constants C and D can be 
determined directly using 



C ^ 



and 



D 



-2Fr{vp) 
(1 - 2,7) 



-2Fi{vp) 
(1 - 27;) 



(13) 



(14) 



These equations are valid if the MPP has been found ex- 
actly, though this will very rarely be the case. Any error 
induced in the magnitude of the Fourier components by 
this effect will be small, but the error induced in the com- 
plex phase (denoted x so as not to be confused with the 
Hamiltonian angle 9) will not be negligible and must be 
corrected. We may do this by observing that in Eq. p2|l 
the constant term and the negative amplitude of the co- 
sine term must be equal. We can define the corrected 
complex angle Xc so that this is the case using 



Xc — arccos 



-F(0) 



2Fr{up) 

such that the corrected Fourier component 

Fc{vp) = |F(t/p)|[cos(xc) + »sin(xc)], 



(15) 



(16) 



is then used in Eqs. and ((Tl|l . 

At this point, in order to keep track of the various 
sine and cosine terms and their uncertainties, we will 
introduce the following notation. When dealing with an 
angle we use = cos(l') and 5A^ to refer to the cosine 
of the angle and its uncertainty respectively. Likewise, 
we define i3$ = sin($) as the sine of the angle giving the 
relationship — ^1 — B\ and A^5A^ — B^5B^. 

As the value of Ok has already been determined, can 
be found from either C or Z?, depending on the value of 



Ok- For instance using 

A^-p 



cos((/) — (3) 
2C/ sin(26'fc) 
C/{AeM 



(17) 



B, 



<t>-0 



sgn{D)Jl-Al_ 



0-/3 



or 



0-/3 



sin((/) — /?) 
i^/sin(0fc) 
D/Be, 



% < 



3tt 



(18) 



A. 



4,-f3 



1 



B2 



depending on the value of Ok, will minimise the effects of 
noise. The angle (/) is then given by 



4> = arccos(A0_^) -t- P, 



(19) 



as expected. 

As the rotation about the axis can only be per- 
formed to the same accuracy as the axis itself is charac- 
terised, there will also be some uncertainty in the angle 
/3. This can be approximated by setting 60^ ~ (5/3, which 
gives the uncertainty 



SA, 



El 

Be^ 



■SA, 



(20) 



in Af 



cos(/3). 



VII. ESTIMATING THE UNCERTAINTY IN 

The uncertainty in cf) will depend on the uncertainty 
in both the original axes characterisation and the noise 
in the Fourier spectrum used to compute C and D. The 
uncertainty in the parameters C and D can be calculated 
using 



and 



sc-" 



2(1 - 27,) 



(1 - 2rj) 



6F^ 



5F' 



2C 



[I - 27^) 



2D 



(1 - 2r;) 



(21) 



(22) 



where SF and Srj are those defined in section^ Here, we 
have ignored the covariance term to simplify the analysis. 
The contribution due to correlated errors is small as the 
calculation of (p depends on three sets of measurements 
(dr, dfe and A^^js) which are independent of each other. 

We can then define the uncertainty in A<^_^ in terms 
of C or as 



"^0-/3 — ^0-/3 



sc_ 
IT 



SAe 

Aa, 



SBg 
Ba, 



(23) 



7 



— /12 



D 



V Be 



Writing the cosine of as 

= cos{(f>) = ApA^^p - BpB^^fi, 
gives the uncertainty relationship 



SAi 



4 



5Al 



5Al^p. 



(24) 



(25) 



(26) 
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VIII. ESTIMATING THE UNCERTAINTY IN 
THE HAMILTONIAN PARAMETERS 

Once the estimates Or, Ok, 4>, oJr and ujk and have been 
found, the Hamiltonians can be estimated using the fol- 
lowing equations, 

Hr = -^(Se^CTj; + Ae^^fTz) 

— Hr^xO'x + Hr^z<^z (27) 

and 

— Hk,xCx + Hk'yCFy + Hk,z<^z, (28) 



where Hj ^ is the i-th component of the j-th Hamiltonian. 
The uncertainty in each of the components of Hr are 
given by 




FIG. 6; An example of the systematic reduction in the un- 
certainty of the Hamiltonian parameters as the number of 
measurements is increased. The error bars are given by three 
times the uncertainty estimate for each point and the solid 
line gives the 'true' value {Hr,x = 0.1, Hr,z ~ 0.05). The es- 
timates are seen to converge to the true value as the number 
of measurements are increased. 



IX. EXAMPLE SIMULATIONS 

To illustrate these ideas and determine the accuracy 
of the parameter estimate and its uncertainty, we sim- 
ulated the measurement procedure on an arbitrary ex- 
ample system, Hr = OAax + O.OSctz. Using an obser- 
vation time tob = 500 and progressively larger numbers 
of measurements, the increase in precision can be ob- 
served. In Fig. El the components Hr.x and Hr^z are 
plotted for increasing numbers of measurements. The er- 
rors bars are given by 3SH which should be equivalent 
to the 3-sigma level and the true value is shown as a 
solid line. As the number of measurements increases, the 
uncertainty reduces and the estimated values converge 
to the true value, as expected. The complete process 
is then simulated using a second example Hamiltonian 
{Hk = 0-Q(Tx + 0.45(Tj, -I- O.ltTz) and similar results are ob- 
tained but with increased uncertainty as the components 
of Hk rely on the measurements of both Hr and Hk , so 
there is more scope for accumulated errors. 

In order to compare the uncertainty calculated using 
the equations in section lVIIII with the expected spread of 
the data, we repeated the simulations of the example sys- 
tem many times with the same number of measurements. 
By looking at the spread of the resulting estimates from 
many experiments and comparing this to the derived un- 
certainty from one experiment we can confirm that the 
uncertainty provides a good bound. Providing a good er- 
ror bound on the Hamiltonian parameters alleviates the 
need to perform characterisation many times to obtain 
good statistics. 
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A. Accuracy of the Uncertainty Estimate 

To measure the distance between the real Hamihonian 
vector d and its estimate d we use the fohowing distance 
metric, 



V ■ 



Id-d| 



(34) 



and a measure of the uncertainties is 



5V. 



Sdl + Sdl + Sdl i^di 



(35) 




0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 



We simulated the characterisation procedure for the 
example system using iob = 500, Ng — 10000 and 
Ne — 50 with a measurement error probability of 10% 
(?7 = 0.1). Fig. EJa) shows a histogram of V for Hr 
over 5000 simulated runs, the average uncertainty ST) 
over 5000 runs is also shown. For this example 98.4% 
of the simulation runs lie within 3(51?, illustrating that 
the uncertainty provides a good bound on the estimated 
parameters. 

The fidelity V for H^. shows a similar distribution, 
though the absolute uncertainty is greater for a given 
number of measurements as more steps are required to 
determine the azimuthal angle (p. Fig. Efb) shows the 
equivalent histogram for determination of the Hamilto- 
nian Hk over 5000 simulated runs. Three times the av- 
erage uncertainty {3ST>) includes 98.7% of the data. The 
intervals for both Vr and Vk are slightly too small as 
a 3-sigma interval should contain approximate 99.7% of 
the data. This discrepancy is due to the effect of cor- 
related errors between the Fourier components SF and 
the uncertainty in the MPP location {Sui). In general, 
as the noise level in the Fourier spectrum increases, the 
width of the peak P{tp) will also increase. This results in 
a small correlation between the uncertainties in d) and 6 
which has not been taken into account. For a given set of 
experimental data, the width of P{tp) and the standard 
deviation of the the noise floor of the Fourier spectrum 
will decrease as the number of measurements increases. 
The relationship between these errors can then be deter- 
mine and will be (in general) non-trivial. The covariance 
can then be calculated, the result of which would be to 
add an additional term to Eqs. I|27|) - (|33|) and therefore 
increasing the overall uncertainty. This additional term 
will be small as the fractional uncertainty in oj is typi- 
cally much smaller than in 9 which implies the covariance 
between them will also be small, relative to the other un- 
certainties. 

The measurement error estimate (?)) is found to be very 
well behaved, with 99.5% of the estimates lying within 
the error bounds, which is very close to what is expected 
for a 3-sigma confidence interval. A histogram of t) is 
shown in Fig. |S1 with 3(5?7 labelled for 5000 runs, each run 
consisting of Nt = 5 x 10^ measurements. 




0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 



FIG. 7: The distribution of D for the estimated (a) Hr and (b) 
Hk over 5000 simulated runs. For these simulations, (a) 98.4% 
and (b) 98.7% of the estimates are found to lie within the 
average uncertainty interval (352?). The absolute uncertainty 
in Hk is greater than for Hr as more steps are required, giving 
a larger accumulated error. 




0.097 0.098 0.099 0.1 0.101 0.102 0.103 



ti (estimated) 

FIG. 8: The distribution of the estimated measurement error 
iff) for 5000 simulated runs. For this simulation, 99.5% of the 
estimates lie within the uncertainty iS^ry. 
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FIG. 9: The average uncertainty SV,. of the estimate for the 
Hamiltonian Hr as a function of total number of measure- 
ments. Each data point is the average of 10 simulation runs. 
The sohd hue shows l/\fN where A'' is the number of mea- 
surements. As the total number of measurements increases, 
the overall precision with which the Hamiltonian is know in- 
creases. For a random measurement error, the achievable pre- 
cision is reduced but still asymptotes to a scaling of one over 
the square-root of the number of measurements. 



B. Scaling behaviour of the Uncertainty 



The usefulness of this technique is ultimately governed 
by how many measurements are required to obtain a 
given precision in the final Hamiltonian estimate. To in- 
vestigate this, the example system was characterised with 
progressively larger numbers of measurements. The aver- 
age of the resulting estimated uncertainty SVr is plotted 
in Fig.|51for several different values of the measurement 
error (77). For increasing numbers of measurements, the 
Hamiltonian estimate gets progressively more accurate, 
as expected. This scaling is approximately proportional 
to 1/VN with the achievable precision reduced by the 
effect of the measurement error. This constant factor 
is effectively a 'penalty' which depends on the measure- 
ment error but is largely independent of the number of 
measurements. 

From this type of analysis we can estimate how many 
measurements are required to achieve a certain precision 
in the final result. Assuming all other factors are negligi- 
ble, the achievable precision scales as one over the square- 
root of the number of measurements. Other factors, such 
as control field fluctuations, will ultimately limit this pro- 
cess. This is easily identified as the achievable precision 
will tend to asymptote to some value which is limited by 
these fluctuations. 



X. IMPLICATIONS FOR SINGLE-QUBIT 
ROTATIONS IN QUANTUM COMPUTING 

In order to be able to perform single-qubit rotations of 
the type required for quantum computing applications, 
a certain level of accuracy is required. The threshold 
theorem for quantum error correction states that if a 
physical error rate of p = 10"'' — 10~^ can be achieved 
then concatenated quantum error correction protocols 
can be implemented successfully for arbitrary precision 
comDutation|12.] . This physical error rate gives the prob- 
ability of a discrete error due to decoherence of the sys- 
tem. The errors introduced due to inaccurate character- 
isation will also contribute, though in a less predictable 
way. Typically, gate operations are assumed to have a 
precision of 10~^ or better but from the previous anal- 
ysis, this would require 10^^ measurements during char- 
acterisation. For a typical measurement readout time of 
1/is this gives an initial characterisation time of approx- 
imately 12 days. 

This turns out to be an overly simplistic view as the 
precision of the gate operations is not equivalent to the 
probability of a discrete error due to decoherence. For 
a single gate rotation around an ideal angle 9, the true 
rotation will be around an angle 6{1 + e) and therefore 
the probability of a discrete error p cx (e)^ where e « 
68/9. Given the previous discussion on the scaling of 69 
with number of characterisation measurements N, the 
probability of discrete error on a single gate operation 
actually scales proportional to N^^, which requires only 
10^ rather than 10^^ measurements. 

As well as errors induced by inaccurate knowledge of 
the Hamiltonian angle (9), errors can also be introduced 
due to an inaccurate rotation frequency or 'over rotation 
error'. In general this will have a similar effect to an 
angle characterisation error as (for small errors) they are 
equivalent. In addition, for the characterisation process 
discussed in this paper, the percentage uncertainty in 
the rotation frequency is typically an order of magnitude 
smaller than the uncertainty in the Hamiltonian angle 
which means that angle errors are the dominant source 
of gate error. 

For multiple gate operations, the probability of a dis- 
crete error scales as np where n is the number of gate 
operation time steps and therefore the number of possi- 
ble error locations [iSj, assuming that errors in different 
qubits are uncorrelated. In the worst case, the rotation 
error accumulates as ne which gives pt = np oa (ne)^, the 
total probability of error for n possible error locations. 
This means its possible (in the worst case) for the uncer- 
tainty in the angle to accumulate over multiple rotations. 
This will not always be the case as certain rotations (such 
as a 2tt rotations) are less susceptible to characterisation 
errors than others and it is possible to get error cancella- 
tion. While this discussion is not new[l2.[T3j. the 1/\/N 
scaling of the achievable precision in 69 highlights the 
very real constraints imposed by the measurement and 
therefore characterisation time of any prospective quan- 
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turn computing proposal. 

Several techniques exist for dealing with characteri- 
sation errors of this kind |3 much of which has 
recently regained interest for QIP applications 0, ^3 ■ 
One such technique, which has been known in the NMR 
literature for some time, is composite pulsing [T^. This 
involves carefully constructing a pulse sequence for a 
given rotation in order to reduce characterisation errors 
in both the angle (off-resonant errors) and the rotation 
frequency (pulse length errors). Recent work by Brown 
et al. [T^ has shown that in-fact, systematic character- 
isation errors can be eliminated to arbitrary order us- 
ing strings of composite pules. For a single imperfect 
gate with fractional error e, the resulting gate error can 
be reduced to ©(e") for arbitrary n using a composite 
pulse sequence whose length scales as n^. Using this or 
similar techniques, we can imagine a trade-off between 
long initial characterisation time (large number of charac- 
terisation measurements) and longer composite pulse se- 
quences for our gate operations (slower operating speed). 
In addition, by choosing fine time sampling (large Ng) 
we can obtain very precise frequency estimates at the ex- 
pense of poor angular resolution due to small numbers 
of ensemble measurements (iVe)- The imprecise angu- 
lar estimate could then be accounted for using composite 
pulsing. Similarly, poor time resolution and large num- 
bers of ensemble measurements will give accurate angle 
estimates at the expense of rotation frequency resolution. 
There may also be situations where it is advantageous to 
precisely characterise some gates and/or qubits but not 
others. 



XI. CONCLUSION 

As the precision and level of complexity of quantum 
control experiments increases, the accuracy to which per- 
tinent system parameters are known must also increase. 
While this is most commonly discussed in the context 
of quantum computing, the ability to precisely measure 
the terms in an arbitrary Hamiltonian has much broader 
application to the study of quantum systems. 

The procedure given here for characterising an arbi- 
trary two-state Hamiltonian has distinct advantages over 
other methods. Given only one measurement axis and as- 
suming the system can be repeatedly initialised in a sin- 
gle state which is not an eigenstate of the Hamiltonian 
to be characterised, the Hamiltonian parameters can be 
determined to arbitrary accuracy. By taking the discrete 
Fourier transform of a series of measurements of the evo- 
lution of the system, the parameters in the Hamiltonian 
can be computed directly from the Fourier components. 

Using signal processing techniques, the uncertainty in 
the Hamiltonian parameters can be estimated and we 
have derived example expressions for these uncertainties. 



If a random measurement error is present, this too can 
be characterised with an uncertainty. This uncertainty 
estimate is found to scale proportionally to one on the 
square-root of the total number of measurements. The 
introduction of measurement error reduces the achievable 
precision by a constant factor which is independent of the 
number of measurements. 

In the laboratory, this procedure can be applied as 
the experiment progresses, giving an increasing more ac- 
curate estimate of the parameters in question. It also 
means that if the response of a Hamiltonian to a given 
input parameter is required, as the input parameter is 
varied, the resulting system can be determine with an 
uncertainty at each point. This enables the usual (non- 
)lincar fitting routines to be applied to the problem to 
find the general response function. 

Being able to accurately characterise a Hamiltonian 
is vitally important if we are to move beyond proof-of- 
concept experiments and build working devices for QIP. 
The trade-off between more accurate initial characterisa- 
tion and more sophisticated gate sequences allows these 
devices to be optimised for a particular application. 
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APPENDIX A: DERIVATION OF THE TIME 
EVOLUTION OF (a^) UNDER AN ARBITRARY 
TWO-STATE HAMILTONIAN 

Given an arbitrary two-state Hamiltonian, we can 
write it in terms of the Pauli matrices using Eq. 
The free evolution of the system under this Hamiltonian 
is given by the operator U (t) = e~*^* which, using a 
generalised de Moivre formula can be rewritten as 



U{t) 



-id„t/2 



I cos 



\d\t_ 

2 



id. (7 sin 



2 



(Al) 

If the system is initially in the state \tp{0)) = \0) {9 = (p = 
0) then (converting to polar coordinates) the evolution of 
the system is given by 
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\i;{t)) = U{t)\m)=e 



_ Adot/2 



COS 



2 



I cos y sin 



2 



|0) + sin ^ sin 



|djt 
2 



(sin^ — icos0)|l) 



(A2) 



r 



The observable in this case is the projection onto the z- 
axis so we will use z = |0)(0| — |1)(1| as the operator 
which gives the expectation value of the z-projection, 

(a,) = (i) = (V'(t)lilV'W)- (A3) 

After cancelling the global phase and rearranging terms, 
this becomes 

(a,) = cos^ f ^) +(cos2 ^-sin^ Q) sin^ . (A4) 



If we set |d| = w, the angular frequency of the precession, 
this gives the time dependence of the z-projection 

z{t) = (cr^) = coswi sin^ 9 + cos^ 6, (A5) 

as expected. 
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